% code used to create Figs. 3 in "An alternative pathway for delivery of 
% power by outer hair cells to cochlear traveling waves", Samaras and Meaud
% plot linear response of the model
% Julien Meaud
% julien.meaud@me.gatech.edu
% ------------------------------------------------------------------------

clear all; close all; clc

r = 2; c = 3;  % number of rows and columns 

setPlottingOptions

panelLabels = {'A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P'};

LWexp=1; % linewdith for experiments
MS = 10; % MarkerSize
FS = 10; % Change FontSize for specific

colV=['b','m','k','c','g','r'];

plotPhaseReBM = 1;

CF = 20e3;  % CF at the plotted location

rEps3 = 1;
modelWithDCs = 1;

LS1 = '-';
LS2 = ':';

TEST = 0;

rDC=logspace(-0.1,3,20);
rRL=logspace(-1,3,20);
nDC = length(rDC);
nRL = length(rRL);

GainAtBP_BM = zeros(nDC,nRL);
powerPCreMax_atMaxNet = zeros(nDC,nRL);
powerDCreMax_atMaxNet = zeros(nDC,nRL);


plotMicromechanics=0; %set to 1 if you want to plot TM and DC

jPlot=0;
for iVar = 1:nDC

    iVar/nDC
    
    tic
    for jVar = 1:nRL
        
        filenamePass = sprintf('ParametricStudy_rDC_rRL_passive/PTLFD_G35_Act1.44eq_rEps3_0_rKdcKohc_%g_Krl_%g_Cdc0_0_25_Cohc0_0_25_Mdc0_001red.mat',rDC(iVar),rRL(jVar));
        filenameAct = sprintf('ParametricStudy_rDC_rRL_active/PTLFD_G35_Act1.85B1_15A_rEps3_1_rDC_%g_Krl_%g_Cdc0_0_25_Cohc0_0_25red.mat',rDC(iVar),rRL(jVar));

        vars_Act = load(filenameAct);
        vars_Pass = load(filenamePass);
       
        freq = vars_Act.StimulusOptionsStruct.freq;
        xMesh = vars_Act.ModelGeneratedVarsStruct.xMesh;
        omega=freq*2*pi;        
        
        vStapes_Act = vars_Act.ModelOutputsStruct.middleEar.vStapes;
        uStapes_Act = vStapes_Act./(2*pi.*freq*1i);
        vStapes_Pas = vars_Pass.ModelOutputsStruct.middleEar.vStapes;
        uStapes_Pas = vStapes_Pas./(2*pi.*freq*1i);
        
        uBM_Act = vars_Act.ModelOutputsStruct.structural.uBM;
        uDC_Act = vars_Act.ModelOutputsStruct.structural.uDC;
        uTMR_Act = vars_Act.ModelOutputsStruct.structural.uTMR;
        uTMB_Act =vars_Act.ModelOutputsStruct.structural.uTMB;
        
        uBM_Pas =vars_Pass.ModelOutputsStruct.structural.uBM;
        uDC_Pas = vars_Pass.ModelOutputsStruct.structural.uDC;
        uTMB_Pas = vars_Pass.ModelOutputsStruct.structural.uTMB;
            
            
        OptionalStruct=vars_Act.OptionalStruct;

        L=vars_Act.ModelGeneratedVarsStruct.xMesh(end);
        
        
        iCFoI = find2(freq,CF);
        [maxAct,iBP] = max(abs(uBM_Act(:,iCFoI)));
        [~,iPeak_pass] = max(abs(uBM_Pas(:,iCFoI)));
        BP = xMesh(iBP);
        
        %% Calculate gain at BP, below BP:
        
        peakBM_act = max(20*log10(abs(uBM_Act(iBP,50:end)./(uStapes_Act(50:end)))));
        peakBM_pass = max(20*log10(abs(uBM_Pas(iBP,50:end)./(uStapes_Pas(50:end)))));

        peakTMB_act = max(20*log10(abs(uTMB_Act(iBP,50:end)./(uStapes_Act(50:end)))));

        
        Q=calculateQ(uBM_Act,freq,iBP,CF,iCFoI,maxAct,6);
        Q6(iVar,jVar)=Q;
      
        Q=calculateQ(uBM_Act,freq,iBP,CF,iCFoI,maxAct,10);
        Q10(iVar,jVar)=Q;
      
        nonlinearityBM_subBF(iVar,jVar)=20*log10(abs(uBM_Act(iBP,int16(iCFoI/2))./uBM_Pas(iBP,int16(iCFoI/2))));
        nonlinearityTMB_subBF(iVar,jVar)=20*log10(abs(uTMB_Act(iBP,int16(iCFoI/2))./uTMB_Pas(iBP,int16(iCFoI/2))));
                
        GainAtBP_BM(iVar,jVar) = peakBM_act - peakBM_pass;

        GainAtBP_TM(iVar,jVar) = peakTMB_act - peakBM_pass;

        maxSigma(iVar,jVar)=max(real(vars_Act.StabilityAnalysisStruct.lambda((~isinf(real(vars_Act.StabilityAnalysisStruct.lambda))))));
        
   
    end
    toc
end

BP

figure(404)
contours =  [0.5 1 1.5 2 2.5 3 3.5 4];
contoursExtend = [0.5 1 1.5 2 2.5 3 3.5 4];

% contours = [-1 2 5 15 25];
% contoursExtend = [-1 2 5 15 25 30];
% Define colormap:
targetMat = contours(1:end-1);%contours(1):1:contours(end)-1;
parulaMat = parula(length(contours)-1);
iContours = find2(targetMat,contoursExtend);
iContours = [iContours, length(targetMat)];
cmapCustom = zeros(length(targetMat),3);
for iC = 1:length(contours)-1
    
    i1 = iContours(iC);
    i2 = iContours(iC+1);
    
    Color = parulaMat(i1,:);
    colorBlock = repmat(Color,i2-i1+1,1);
    
    cmapCustom(i1:i2,:) = colorBlock;
end

subplot(1,2,1)
[C,h]=contourf(rDC,rRL,transpose(Q6),contours); % transpose needed if Kdc is x-axis, Kohc is y-axis
h.LineStyle='none';

set(gca,'xscale','log')
set(gca,'yscale','log')
xlabel('r_{dc}')
ylabel('r_{rl}')
% clabel(C,h)
cbar = contourcbar;
cbar.Ticks = contours(2:end);
colormap(cmapCustom)
caxis([contours(1) contours(end)]);%[-1 42])
cbar.Label.String = 'Q_{6}';
cbar.Location = 'northoutside';
axis tight


figure(4)
subplot(r,c,2*c)
[C,h]=contourf(rDC,rRL,transpose(Q10),contours); % transpose needed if Kdc is x-axis, Kohc is y-axis
h.LineStyle='none';
hold on

set(gca,'xscale','log')
set(gca,'yscale','log')
xlabel('r_{dc}')
ylabel('r_{rl}')
% clabel(C,h)
cbar = contourcbar;
cbar.Ticks = contours(2:end);
colormap(cmapCustom)
caxis([contours(1) contours(end)]);%[-1 42])

cbar.Location = 'eastoutside';
axis tight
text(10,3000,'Q_{10}','FontSize',12,'FontWeight','Bold')


if max(max(maxSigma))>0
    [CStab,hStab]=contour(rDC,rRL,maxSigma.',[0,1e8]);
    hStab.LineWidth = 2;
    hStab.LineColor='k';
    
    Ni=CStab(2,1);
    xData=CStab(1,2:Ni+1);
    yData=CStab(2,2:Ni+1);
    
    %patch([x fliplr(x)], [G1 min(ylim)*ones(size(G1))], 'r')        % Below Lower Curve
    patch([xData fliplr(xData)], [yData max(ylim)*ones(size(yData))], 'w','EdgeColor','none')        % Above Upper Curve
    alpha(0.5)

    text(25,200,'Unstable','Fontsize',10,'Color','k','Fontweight','bold')
    text(5,2,'Stable','Fontsize',10,'Color','k','Fontweight','bold')

end
%hold off

text(1.2,500,'F','Fontsize',16,'Color','w','Fontweight','bold')
text(5,25,'M4','Fontsize',10,'Color','k','Fontweight','bold')
text(180,12,'M3','Fontsize',10,'Color','k','Fontweight','bold')
text(6,0.36,'M2','Fontsize',10,'Color','k','Fontweight','bold')
text(180,0.2,'M1','Fontsize',10,'Color','k','Fontweight','bold')
        
%% Plot surface plot of gain  --------------------------------------------

contours =  [-4 1 7 15 22 30 38 46];
contoursExtend = [-4 1 7 15 22 30 38 46];

% contours = [-1 2 5 15 25];
% contoursExtend = [-1 2 5 15 25 30];
% Define colormap:
targetMat = contours(1):1:contours(end)-1;
parulaMat = parula(length(contours(1):1:contours(end))-1);
iContours = find2(targetMat,contoursExtend);
iContours = [iContours, length(targetMat)];
cmapCustom = zeros(length(targetMat),3);

for iC = 1:length(contours)-1
    
    i1 = iContours(iC);
    i2 = iContours(iC+1);
    
    Color = parulaMat(i1,:);
    colorBlock = repmat(Color,i2-i1+1,1);
    
    cmapCustom(i1:i2,:) = colorBlock;
end
xticks([0.1 1 10 100 1000])
yticks([0.1 1 10 100 1000])

figure(4)
subplot(r,c,c)
[C,h]=contourf(rDC,rRL,GainAtBP_BM.',contours); % transpose needed if Kdc is x-axis, Kohc is y-axis
h.LineStyle='none';
set(gca,'xscale','log')
set(gca,'yscale','log')
xlabel('r_{dc}')
ylabel('r_{rl}')
% clabel(C,h)
cbar = contourcbar;
cbar.Ticks = contours(2:end);
colormap(cmapCustom)
caxis([contours(1) contours(end)]);%[-1 42])
%cbar.Label.String = 'G^{peak}_{act/pass} (dB)';
cbar.Location = 'eastoutside';
axis tight

shading interp
set(gca,'Fontsize',10)
ylim([rRL(1) rRL(end)])
xlim([rDC(1) rDC(end)])

xticks([0.1 1 10 100 1000])
yticks([0.1 1 10 100 1000])

text(3,4000,'G^{peak}_{act/pass} (dB)','FontSize',12,'FontWeight','Bold')

[~,iMaxGainVsKdc] = max(GainAtBP_BM.');
hold on

% add the stability limit contour -----------------------------------------
if max(max(maxSigma))>0
    [CStab,hStab]=contour(rDC,rRL,maxSigma.',[0,1e8]);
    hStab.LineWidth = 2;
    hStab.LineColor='k';
    
    Ni=CStab(2,1);
    xData=CStab(1,2:Ni+1);
    yData=CStab(2,2:Ni+1);
    
    %patch([x fliplr(x)], [G1 min(ylim)*ones(size(G1))], 'r')        % Below Lower Curve
    patch([xData fliplr(xData)], [yData max(ylim)*ones(size(yData))], 'w','EdgeColor','none')        % Above Upper Curve
    alpha(0.5)

    text(25,200,'Unstable','Fontsize',10,'Color','k','Fontweight','bold')
    text(5,2,'Stable','Fontsize',10,'Color','k','Fontweight','bold')
end
%hold off

text(1.2,500,'E','Fontsize',16,'Color','w','Fontweight','bold')
text(5,25,'M4','Fontsize',10,'Color','k','Fontweight','bold')
text(180,12,'M3','Fontsize',10,'Color','k','Fontweight','bold')
text(6,0.36,'M2','Fontsize',10,'Color','k','Fontweight','bold')
text(180,0.2,'M1','Fontsize',10,'Color','k','Fontweight','bold')



figure(402)
[C,h]=contourf(rDC,rRL,GainAtBP_TM.',contours); % transpose needed if Kdc is x-axis, Kohc is y-axis
h.LineStyle='none';
set(gca,'xscale','log')
set(gca,'yscale','log')
xlabel('r_{dc}')
ylabel('r_{rl}')
% clabel(C,h)
cbar = contourcbar;
cbar.Ticks = contours(2:end);
colormap(cmapCustom)
caxis([contours(1) contours(end)]);%[-1 42])
cbar.Label.String = 'TM G^{peak}_{act/pass} (dB)';
cbar.Location = 'northoutside';
axis tight

shading interp
set(gca,'Fontsize',10)
ylim([rRL(1) rRL(end)])
xlim([rDC(1) rDC(end)])

xticks([0.1 1 10 100 1000])
yticks([0.1 1 10 100 1000])


[~,iMaxGainVsKdc] = max(GainAtBP_BM.');
hold on

% add the stability limit contour -----------------------------------------
if max(max(maxSigma))>0
    [CStab,hStab]=contour(rDC,rRL,maxSigma.',[0,1e8]);
    hStab.LineWidth = 2;
    hStab.LineColor='k';
    
    Ni=CStab(2,1);
    xData=CStab(1,2:Ni+1);
    yData=CStab(2,2:Ni+1);
    
    %patch([x fliplr(x)], [G1 min(ylim)*ones(size(G1))], 'r')        % Below Lower Curve
    patch([xData fliplr(xData)], [yData max(ylim)*ones(size(yData))], 'w','EdgeColor','none')        % Above Upper Curve
    alpha(0.5)

    text(25,200,'Unstable','Fontsize',10,'Color','k','Fontweight','bold')
    text(5,2,'Stable','Fontsize',10,'Color','k','Fontweight','bold')

end







% surface plot of the maximum real part of the eigenvalues ----------------
figure(401)
surf(rDC,rRL,maxSigma.','EdgeColor','none'); 
hold on
[C,h]=contour(rDC,rRL,maxSigma.',[0,1e8]);
h.LineWidth = 2;
h.LineColor='k';
view(0,90)
set(gca,'xscale','log')
set(gca,'yscale','log')
xlabel('r_{dc}')
ylabel('r_{rl}')
cbar = colorbar;

cbar.Label.String = 'max \sigma';
cbar.Location = 'northoutside';
axis tight

set(gca,'Fontsize',10)
ylim([rRL(1) rRL(end)])
xlim([rDC(1) rDC(end)])

xticks([0.1 1 10 100 1000])
yticks([0.1 1 10 100 1000])


%% Plot BM responses for models M1, M2, M3, M4 and M5  --------------------
%figure(2)
LW = 2;
LS1 = '-';
LS2 = ':';

titleCell = {'r_{dc}=1, r_{rl}=12.3', 'r_{dc}=10, r_{rl}=12.3', 'r_{dc}=1000, r_{rl}=12.3',...
    'r_{dc}=1, r_{rl}=0.14', 'r_{dc}=10, r_{rl}=0.14', 'r_{dc}=1000, r_{rl}=0.14'};

    filenamesAct_BL_Cell ={'PTLFD_G35_Act1.85B1_15A_rEps3_1_rKdcKohc_1000_Krl_BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat',...                           
    'PTLFD_G35_Act1.85B1_15A_rEps3_1_rKdcKohc_1000_Krl_90BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat',...  
    'PTLFD_G35_Act1.85B1_15A_rEps3_1_rKdcKohc_5_Krl_BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat',... 
    'PTLFD_G35_Act1.85B1_15A_rEps3_1_rKdcKohc_5_Krl_90BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat',...
    'PTLFD_G35_Act1.85B1_15A_rEps3_1_rKdcKohc_5_Krl_365BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat'};
  
     filenamesAct_Cell = {'PTLFD_G35_Act1.85B1_15A_rEps3_1_rKdcKohc_1000_Krl_BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat',...
     'PTLFD_G35_Act1.65B1_15A_rEps3_1_rKdcKohc_1000_Krl_90BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat',...  
    'PTLFD_G35_Act3.15B1_15A_rEps3_1_rKdcKohc_5_Krl_BL_Cdc0_0_25_Cohc0_0_25.mat',... 
    'PTLFD_G35_Act3.0B1_15A_rEps3_1_rKdcKohc_5_Krl_90BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat',...
    'PTLFD_G35_Act3.78B1_15A_rEps3_1_rKdcKohc_5_Krl_365BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat'};

filenamesPass_Cell = {'PTLFD_G35_Act1.44eq_rEps3_0_rKdcKohc_1000_Krl_BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat',...
     'PTLFD_G35_Act1.44eq_rEps3_0_rKdcKohc_1000_Krl_90BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat',...  
    'PTLFD_G35_Act1.44eq_rEps3_0_rKdcKohc_5_Krl_BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat',... 
    'PTLFD_G35_Act1.44eq_rEps3_0_rKdcKohc_5_Krl_90BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat',...
    'PTLFD_G35_Act1.85B1_15A_rEps3_0_rKdcKohc_5_Krl_365BL_Cdc0_0_25_Cohc0_0_25_Mdc0_001.mat'};


%cd Models_M1M2M3M4_FreqDomain

iPlotVec = [1 4 2 5 6];
iPlotVecPanoramic = [1 3 2 4 6];
modelTitle=['1','2','0','3','4','5'];
iPanelLabels = [4 3 2 1 5];
iFilesVec = [1 2 3 4 5];
titleStr=['A','B','C','C','D','E','F'];


for iPlot = [1,4,3,2]  % add 5 if want to plot M5
    filenameAct_BL = ['Models_M1M2M3M4_FreqDomain/',filenamesAct_BL_Cell{iPlot}];
    filenameAct = ['Models_M1M2M3M4_FreqDomain/',filenamesAct_Cell{iPlot}];
    filenamePass = ['Models_M1M2M3M4_FreqDomain/',filenamesPass_Cell{iPlot}];    
    
    vars_Act_BL= load(filenameAct_BL);
    vars_Act= load(filenameAct);
    vars_Pass = load(filenamePass);
    
    freq = vars_Act.StimulusOptionsStruct.freq;
    xMesh = vars_Act.ModelGeneratedVarsStruct.xMesh;

    lambda_Act=vars_Act.StabilityAnalysisStruct.lambda;
    lambda_Act_BL=vars_Act_BL.StabilityAnalysisStruct.lambda;
    
    vStapes_Act_BL = vars_Act_BL.ModelOutputsStruct.middleEar.vStapes;
    uStapes_Act_BL = vStapes_Act_BL./(2*pi.*freq*1i);
    vStapes_Act = vars_Act.ModelOutputsStruct.middleEar.vStapes;
    uStapes_Act = vStapes_Act./(2*pi.*freq*1i);
    vStapes_Pas = vars_Pass.ModelOutputsStruct.middleEar.vStapes;
    uStapes_Pas = vStapes_Pas./(2*pi.*freq*1i);

    actLevel=vars_Act.OptionalStruct.GaMax(1)./vars_Act_BL.OptionalStruct.GaMax(1);   

    rDC=vars_Act_BL.OptionalStruct.calibration.cofKdc/vars_Act_BL.OptionalStruct.calibration.cofKohc;
    rRL=vars_Act_BL.OptionalStruct.calibration.cofKrl/vars_Act_BL.OptionalStruct.calibration.cofKohc;

    figure(4)
    subplot(r,c,c);
    plot(rDC,rRL,'x','MarkerSize',8,'Linewidth',2,'Color','k') 

    subplot(r,c,2*c);
    plot(rDC,rRL,'x','MarkerSize',8,'Linewidth',2,'Color','k') 

    uBM_Act_BL = vars_Act_BL.ModelOutputsStruct.structural.uBM;

    uBM_Act = vars_Act.ModelOutputsStruct.structural.uBM;
    uDC_Act = vars_Act.ModelOutputsStruct.structural.uDC;
    uTMB_Act = vars_Act.ModelOutputsStruct.structural.uTMB;   

    uBM_Pass = vars_Pass.ModelOutputsStruct.structural.uBM;
    uDC_Pass = vars_Pass.ModelOutputsStruct.structural.uDC;
    uTMB_Pass = vars_Pass.ModelOutputsStruct.structural.uTMB;

    nonlinearityBMSubBF(iPlot)=20*log10(abs(uBM_Act(iBP,int16(iCFoI/2))./uBM_Pas(iBP,int16(iCFoI/2))));
    nonlinearityTMSubBF(iPlot)=20*log10(abs(uTMB_Act(iBP,int16(iCFoI/2))./uTMB_Pas(iBP,int16(iCFoI/2))));
    nonlinearityDCSubBF(iPlot)=20*log10(abs(uDC_Act(iBP,int16(iCFoI/2))./uDC_Pas(iBP,int16(iCFoI/2))));
   
    
    stimLevel =  vars_Act.StimulusOptionsStruct.stimulusAmplitude.level;
    stimulusPa = convertSPLtoPressureInPascals(stimLevel);
    
  
 
        iCFoI = find2(freq,CF);

        if CF==20000            
            [maxAct,iBP] = max(abs(uBM_Act(:,iCFoI)));
            [maxActBL,iBP_BL] = max(abs(uBM_Act_BL(:,iCFoI)));
            [~,iPeak_pass] = max(abs(uBM_Pass(:,iCFoI)));

            [Q_BL,freq1_BL,freq2_BL]=calculateQ(uBM_Act_BL,freq,iBP_BL,CF,iCFoI,maxActBL,6);
            Qmodel_BL(iPlot)=Q_BL;
    
            [Q,freq1,freq2]=calculateQ(uBM_Act,freq,iBP,CF,iCFoI,maxAct,6);
            Qmodel(iPlot)=Q;

            GainAtBP_BM = 20*log10(abs(uBM_Act(iBP,iCFoI))) - 20*log10(abs(uBM_Pass(iPeak_pass,iCFoI)));
            GainAtBP_BL_BM = 20*log10(abs(uBM_Act_BL(iBP_BL,iCFoI))) - 20*log10(abs(uBM_Pass(iPeak_pass,iCFoI)))  ;

            trueCF=20000;

        elseif CF==5000

            for iX=1:length(xMesh)
                [~,iCF]=max(abs(uBM_Act(iX,50:end)));
                CFmap(iX)=freq(iCF+50-1);
            end


            iBP=find2(CFmap,CF);
            
            
            for iX=1:length(xMesh)
                [~,iCF]=max(abs(uBM_Act_BL(iX,50:end)));
                CFmap_BL(iX)=freq(iCF+50-1);
            end

            
            iBP_BL=find2(CFmap_BL,CF);

        end
    
    
        BP(iPlot)=xMesh(iBP);
        
        figure(4)
        subplot(r,c,iPlotVec(iPlot))
        set(gca,'fontsize',10)        
        plot(freq./1e3,20*log10(abs(uBM_Act_BL(iBP,:))./max(abs(uBM_Pass(iBP,100:end)))),'Color','b','LineWidth',LW-0.5,'Linestyle',LS1)
        hold on
         if iPlot>0
            plot(freq./1e3,20*log10(abs(uBM_Act(iBP,:))./max(abs(uBM_Pass(iBP,100:end)))),'Color','b','LineWidth',LW+0.5,'Linestyle',LS1)        
            text(4*CF/20000,33,[num2str(actLevel,'%.2f'),'G_{hb,BL}^{max}'],'Fontsize',9) 
            hold on
        
            % vline([freq1/1e3 freq1/1e3])
            % vline([freq2/1e3 freq2/1e3])
        
            if plotMicromechanics==1
                plot(freq./1e3,20*log10(abs(uDC_Act(iBP,:))./max(abs(uBM_Pass(iBP,100:end)))),'Color','m','LineWidth',LW+0.5,'Linestyle',LS1) % BM conversion from cm to nm (1e7) and normalized by stimulus in pascals
                plot(freq./1e3,20*log10(abs(uTMB_Act(iBP,:))./max(abs(uBM_Pass(iBP,100:end)))),'Color','g','LineWidth',LW+0.5,'Linestyle',LS1) % BM conversion from cm to nm (1e7) and normalized by stimulus in pascals
            end
        end
        plot(freq./1e3,20*log10(abs(uBM_Pass(iBP,:))./max(abs(uBM_Pass(iBP,100:end)))),'Color','b','LineWidth',LW+0.5,'Linestyle',LS2)
         if plotMicromechanics==1
             plot(freq./1e3,20*log10(abs(uDC_Pass(iBP,:))./max(abs(uBM_Pass(iBP,100:end)))),'Color','m','LineWidth',LW+0.5,'Linestyle',LS2)
            plot(freq./1e3,20*log10(abs(uTMB_Pass(iBP,:))./max(abs(uBM_Pass(iBP,100:end)))),'Color','g','LineWidth',LW+0.5,'Linestyle',LS2)
         end
        
         xlim([3 30]*CF/20000)
         ylim([-50 45])
      
   
    vline(CF/1000,'k--')
    
    % addd ticks, tick labels and xlabel ----------------
    if CF==20000
        xticks([5:5:30])
    elseif CF==5000
        xticks([1:1:10])
    end

    if any(iPlotVec(iPlot)==[4 5])
        if CF==20000
            xticklabels([5:5:30])
        elseif CF==5000
            xticklabels([1:1:10])
        end        
        xlabel('Freq. (kHz)','Fontsize',12)
    else
        xticklabels([]);
    end
    yticks([-40 -20 0 20 40])
    if any(iPlotVec(iPlot)==[1 4])
       yticklabels([-40,-20,0,20,40])
    else
     yticklabels([])
    end
    
    text(15*CF/20000,55,['M',modelTitle(iPlotVec(iPlot))],'Fontsize',12)
    text(4*CF/20000,-40,titleStr(iPlotVec(iPlot)),'Fontsize',16,'FontWeight','bold')

    figure(40)
    subplot(2,5,iPlot)
     for iEig=1:length(lambda_Act_BL)
        if real(lambda_Act_BL(iEig))>0
            col='r';
        else
            col='b';
        end

        undampingRatio=real(lambda_Act_BL(iEig))./sqrt(imag(lambda_Act_BL(iEig)).^2+real(lambda_Act_BL(iEig)).^2);
        plot(imag(lambda_Act_BL(iEig))/1000,undampingRatio,'LineStyle','none','Marker','.','Color',col)
        hold on
     end    
   xlim([0.5 100])
    xlabel('Frequency (kHz)')
    ylabel('\sigma (ms^{-1})')
    ylim([-1 0.1])
    line([0.5 100],[0 0],'Color','k','LineStyle','--')
    title(['M',int2str(iPlot),', baseline MET'])
   
    if iPlot>0
        subplot(2,5,5+iPlot)
        for iEig=1:length(lambda_Act)
            if real(lambda_Act(iEig))>0
                col='r';
            else
                col='b';
            end
    
            undampingRatio=real(lambda_Act(iEig))./sqrt(imag(lambda_Act(iEig)).^2+real(lambda_Act(iEig)).^2);
            plot(imag(lambda_Act(iEig))/1000,undampingRatio,'LineStyle','none','Marker','.','Color',col)
            hold on
                
        end
        xlim([0.5 100])
        xlabel('Frequency (kHz)')
        ylabel('\sigma (ms^{-1})')
        ylim([-1 0.1])
        line([0.1 100],[0 0],'Color','k','LineStyle','--')
        title(['M',int2str(iPlot),', adjusted MET'])
    end

  % plot phase of BM displacement -----------------------------------------
    figure(41)
    subplot(r,c,iPlotVec(iPlot))
    set(gca,'fontsize',10) 
    if any(iPlot==[2:6])
        plot(freq./1e3,unwrap(angle(uBM_Act(iBP,:)))/(2*pi),'Color','b','LineWidth',LW-0.5,'Linestyle',LS1) % BM conversion from cm to nm (1e7) and normalized by stimulus in pascals
        text(4*CF/20000,33,[num2str(actLevel,'%.2f'),'G_{hb,BL}^{max}'],'Fontsize',9) 
    end
    hold on
    plot(freq./1e3,unwrap(angle(uBM_Act_BL(iBP,:)))/(2*pi),'Color','b','LineWidth',LW+0.5,'Linestyle',LS1)
    plot(freq./1e3,unwrap(angle(uBM_Pass(iBP,:)))/(2*pi),'Color','b','LineWidth',LW,'Linestyle',LS2)    
    xlim([3 30])
    ylim([-5 1])
    ylabel('BM phase (cycles)')
    
    % plot wavelength traveling wave on BM -----------------------------------------
    figure(42)
    subplot(r,c,iPlotVec(iPlot))
    set(gca,'fontsize',10) 
    if any(iPlot==[2:6])
        for iFreq=1:length(freq)            
            wavelength=-2*pi./(diff(unwrap(angle(uBM_Act(:,iFreq))))./transpose(diff(xMesh)));
            wavelengthM_Act(:,iFreq)=[wavelength(1);wavelength];            
        end
        semilogy(freq./1e3, wavelengthM_Act(iBP,:)*1e4,'Color','b','LineWidth',LW-0.5,'Linestyle',LS1) 
        % text(7,33,[num2str(actLevel,'%.2f'),'G_{hb,BL}^{max}'],'Fontsize',9) 
    end
    hold on
    for iFreq=1:length(freq)           
            wavelength=-2*pi./(diff(unwrap(angle(uBM_Act_BL(:,iFreq))))./transpose(diff(xMesh)));
            wavelengthM_Act_BL(:,iFreq)=[wavelength(1);wavelength];
   
        end
     semilogy(freq./1e3, wavelengthM_Act_BL(iBP,:)*1e4,'Color','b','LineWidth',LW+0.5,'Linestyle',LS1) 

    for iFreq=1:length(freq)           
        wavelength=-2*pi./(diff(unwrap(angle(uBM_Pass(:,iFreq))))./transpose(diff(xMesh)));
        wavelengthM_Pass(:,iFreq)=[wavelength(1);wavelength];
    
    end
     semilogy(freq./1e3, wavelengthM_Pass(iBP,:)*1e4,'Color','b','LineWidth',LW,'Linestyle',LS2) 
    
    xlim([3 30])
    ylabel('Wavelength ({\mu}m)')
    set(gca,'YScale','log')
     
  % plot response  at multiple longitudinal locations
    figure(60)
    xVec=[41 81 121 161 201 241 281 321 ];    
   transp = linspace(0.3,1,length(xVec));
    for iX=1:length(xVec)
         [~,indexBF]=max(abs(uBM_Act(xVec(iX),51:end)));
        BF=freq(indexBF+51);

         subplot(3,6,iPlotVecPanoramic(iPlot))
         plot(freq/BF,20*log10(abs(uBM_Act(xVec(iX),:)./max(abs(uBM_Pass(xVec(iX),:))))),'Color',[blue_CB transp(iX)],'LineWidth',2)
         hold on
         plot(freq/BF,20*log10(abs(uBM_Pass(xVec(iX),:)./max(abs(uBM_Pass(xVec(iX),:))))),'Color',[blue_CB transp(iX)],'LineStyle','--','LineWidth',2)
         xlabel('Norm. frequency')
         xlim([0.3 1.5]);ylim([-20 50])
         ylabel('BM disp. re max pass (dB)')  
    
         subplot(3,6,6+iPlotVecPanoramic(iPlot))
         xlabel('Norm. frequency')
         xlim([0.3 1.5]);ylim([-20 50])
         ylabel('DC disp. re BM disp. (dB)')

         subplot(3,6,12+iPlotVecPanoramic(iPlot))
         xlabel('Norm. frequency ')
         xlim([0.3 1.5]);ylim([-0.5 0.5])
         ylabel('Phase DC re BM (cycle)')
    end

     for iX=1:length(xMesh)
            [maxUbmAct,indexBF]=max(abs(uBM_Act(iX,50:end)));
            BF(iX)=freq(49+indexBF);
            UbmAct_peak(iX)=maxUbmAct;

            [maxUbmAct_BL,indexBF_BL]=max(abs(uBM_Act_BL(iX,50:end)));
            BF_BL(iX)=freq(49+indexBF_BL);
            UbmAct_BL_peak(iX)=maxUbmAct_BL;

            [maxUbmPass,indexBFpass]=max(abs(uBM_Pass(iX,50:end)));
            BFpass(iX)=freq(49+indexBFpass);
            UbmPass_peak(iX)=maxUbmPass;

            activeToPassiveGain(iX)=20*log10(abs(maxUbmAct./maxUbmPass));
            activeToPassiveGain_BL(iX)=20*log10(abs(maxUbmAct_BL./maxUbmPass));
     end
     figure(62)
     subplot(2,2,1)
     plot(xMesh,activeToPassiveGain,'Color',colV(iPlot),'DisplayName',['Adjusted, M',modelTitle(iPlotVec(iPlot))])
     hold on
     plot(xMesh,activeToPassiveGain_BL,'Color',colV(iPlot),'LineStyle','--','DisplayName',['BL, M',modelTitle(iPlotVec(iPlot))])
     xlabel('x (cm)');ylabel('G_{act/pass}^{peak} (dB)')
     legend

     subplot(2,2,2)
     semilogx(BF/1000,activeToPassiveGain,'Color',colV(iPlot))
     hold on
      semilogx(BF_BL/1000,activeToPassiveGain_BL,'Color',colV(iPlot),'LineStyle','--')
     xlabel('BF (kHz)');ylabel('G_{act/pass}^{peak} (dB)')
     xlim([0.5 80])

    subplot(2,2,3)
    semilogy(xMesh,BF/1000,'Color',colV(iPlot))
    hold on
    semilogy(xMesh,BF_BL/1000,'Color',colV(iPlot),'LineStyle','--')
    semilogy(xMesh,BFpass/1000,'Color',colV(iPlot),'LineStyle',':')
    xlabel('x (cm)');ylabel('BF')
    ylim([0.5 80])
   % xlim([0.50 80])

end

fig4=figure(4);
subplot(r,c,1)
text(-4*CF/20000,-165,'BM displacement re passive peak (dB)','rotation',90,'Fontsize',12)
subplot(r,c,2)
legend(['Act. w/',char(10),'baseline MET'],['Act. w/',char(10),'adjusted MET'],'Passive','FontSize',10,'Orientation','horizontal')
legend boxoff

if r==2
  %  set(subplot(r,c,3),'Position',[0.73 0.3253 0.25 0.3991])
    
  set(subplot(r,c,3),'Position',[0.73 0.6 0.2 0.3]);set(gca,'Fontsize',10)
  set(subplot(r,c,6),'Position',[0.73 0.1 0.2 0.3]);set(gca,'Fontsize',10)

    set(subplot(r,c,1),'Position',[0.08 0.65 0.25 0.3])
    set(subplot(r,c,2),'Position',[0.38 0.65 0.25 0.3])
    
    set(subplot(r,c,4),'Position',[0.08 0.25 0.25 0.3])
    set(subplot(r,c,5),'Position',[0.38 0.25 0.25 0.3])

    set(fig4,'units','inches','position',[2 2 6 4])

end

